!  Given the Arnoldi factorization
!     A*V_{k} - V_{k}*H_{k} = r_{k+p}*e_{k+p}^T,
!  apply NP implicit shifts resulting in
!     A*(V_{k}*Q) - (V_{k}*Q)*(Q^T* H_{k}*Q) = r_{k+p}*e_{k+p}^T * Q
!  where Q is an orthogonal matrix which is the product of rotations
!  and reflections resulting from the NP bulge change sweeps.
!  The updated Arnoldi factorization becomes:
!
!     A*VNEW_{k} - VNEW_{k}*HNEW_{k} = rnew_{k}*e_{k}^T.
!\Arguments
!  N       Integer.  (INPUT)
!          Problem size, i.e. size of matrix A.
!
!  KEV     Integer.  (INPUT/OUTPUT)
!          KEV+NP is the size of the input matrix H.
!          KEV is the size of the updated matrix HNEW.
!
!  NP      Integer.  (INPUT)
!          Number of implicit shifts to be applied.
!
!  SHIFT   complex(8) array of length NP.  (INPUT)
!          The shifts to be applied.
!
!  V       complex(8) N by (KEV+NP) array.  (INPUT/OUTPUT)
!          On INPUT, V contains the current KEV+NP Arnoldi vectors.
!          On OUTPUT, V contains the updated KEV Arnoldi vectors
!          in the first KEV columns of V.
!
!  LDV     Integer.  (INPUT)
!          Leading dimension of V exactly as declared in the calling
!          program.
!
!  H       complex(8) (KEV+NP) by (KEV+NP) array.  (INPUT/OUTPUT)
!          On INPUT, H contains the current KEV+NP by KEV+NP upper
!          Hessenberg matrix of the Arnoldi factorization.
!          On OUTPUT, H contains the updated KEV by KEV upper Hessenberg
!          matrix in the KEV leading submatrix.
!
!  LDH     Integer.  (INPUT)
!          Leading dimension of H exactly as declared in the calling
!          program.
!
!  RESID   complex(8) array of length N.  (INPUT/OUTPUT)
!          On INPUT, RESID contains the the residual vector r_{k+p}.
!          On OUTPUT, RESID is the update residual vector rnew_{k}
!          in the first KEV locations.
!
!  Q       complex(8) KEV+NP by KEV+NP work array.  (WORKSPACE)
!          Work array used to accumulate the rotations and reflections
!          during the bulge chase sweep.
!
!  LDQ     Integer.  (INPUT)
!          Leading dimension of Q exactly as declared in the calling
!          program.
!
!  WORKL   complex(8) work array of length (KEV+NP).  (WORKSPACE)
!          Private (replicated) array on each PE or array allocated on
!          the front end.
!
!  WORKD   complex(8) work array of length 2*N.  (WORKSPACE)
!          Distributed array used in the application of the accumulated
!          orthogonal matrix Q.
!-----------------------------------------------------------------------
subroutine cnapps(n, kev, np, shift, v, ldv, h, ncv, resid, q, workd)
    implicit none
    integer kev, ncv, ldv, n, np
    complex h(ncv, kev + np), resid(n), shift(np), v(ldv, kev + np), &
            q(ncv, kev + np), workd(2 * n)
    ! local variables
    complex, parameter :: one = (1.0, 0.0), zero = (0.0, 0.0)
    integer    i, iend, istart, j, jj, kplusp
    logical    first
    complex cdum, f, g, h11, h21, r, s, sigma, t, workl(kev + np)
    real    c, ovfl, smlnum, ulp, unfl, tst1, re, im
    save       first, ovfl, smlnum, ulp, unfl
    ! External Functions
    real, external :: clanhs, slamch, slapy2
    ! Statement Functions
    real   cabs1
    cabs1(cdum) = abs(real(cdum)) + abs(aimag(cdum))
    ! Data statments
    data       first / .true. /

    ! Executable Statements
    if (first) then
        ! Set machine-dependent constants for the stopping criterion.
        ! If norm(H) <= sqrt(OVFL), overflow should not occur.
        ! REFERENCE: LAPACK subroutine zlahqr
        unfl = slamch('safe minimum')
        ovfl = real(1.0 / unfl)
        call slabad(unfl, ovfl)
        ulp = slamch('precision')
        smlnum = unfl * (n / ulp)
        first = .false.
    end if
    kplusp = kev + np
    ! Initialize Q to the identity to accumulate the rotations and reflections
    call claset ('All', kplusp, kplusp, zero, one, q, ncv)
    ! Quick return if there are no shifts to apply
    if (np == 0) return
    ! Chase the bulge with the application of each implicit shift.
    ! Each shift is applied to the whole matrix including each block.
    do jj = 1, np
        sigma = shift(jj)
        istart = 1
20      continue
        do i = istart, kplusp - 1
            ! Check for splitting and deflation. Use a standard test as in the QR algorithm
            ! REFERENCE: LAPACK subroutine zlahqr
            tst1 = cabs1(h(i, i)) + cabs1(h(i + 1, i + 1))
            if(tst1 == 0.0) tst1 = clanhs('1', kplusp - jj + 1, h, ncv, workl)
            if (abs(real(h(i + 1, i))) <= max(ulp * tst1, smlnum))  then
                iend = i
                h(i + 1, i) = 0.0
                go to 40
            end if
        end do
        iend = kplusp
40      continue

        ! No reason to apply a shift to block of order 1 or if the current
        ! block starts after the point of compression since we'll discard this stuff
        if (istart == iend .or. istart > kev) go to 100
        h11 = h(istart, istart)
        h21 = h(istart + 1, istart)
        f = h11 - sigma
        g = h21

        do i = istart, iend - 1
            ! Construct the plane rotation G to zero out the bulge
            call clartg (f, g, c, s, r)
            if (i > istart) then
                h(i, i - 1) = r
                h(i + 1, i - 1) = 0.0
            end if
            ! Apply rotation to the left of H;  H <- G'*H
            do j = i, kplusp
                t = c * h(i, j) + s * h(i + 1, j)
                h(i + 1, j) = -conjg(s) * h(i, j) + c * h(i + 1, j)
                h(i, j) = t
            end do
            ! Apply rotation to the right of H;  H <- H*G
            do j = 1, min(i + 2, iend)
                t = c * h(j, i) + conjg(s) * h(j, i + 1)
                h(j, i + 1) = -s * h(j, i) + c * h(j, i + 1)
                h(j, i) = t
            end do
            ! Accumulate the rotation in the matrix Q;  Q <- Q*G'
            do j = 1, min(i + jj, kplusp)
                t = c * q(j, i) + conjg(s) * q(j, i + 1)
                q(j, i + 1) = - s * q(j, i) + c * q(j, i + 1)
                q(j, i) = t
            end do
            ! Prepare for next rotation
            if (i < iend - 1) then
                f = h(i + 1, i)
                g = h(i + 2, i)
            end if
        end do
        ! Finished applying the shift.
100     continue
        ! Apply the same shift to the next block if there is any.
        istart = iend + 1
        if (iend < kplusp) go to 20
        ! Loop back to the top to get the next shift.
    end do
    ! Perform a similarity transformation that makes sure that
    ! the compressed H will have non-negative real subdiagonal elements.
    do j = 1, kev
        re = real(h(j + 1, j))
	im = aimag(h(j + 1, j))
        if (re < 0.0 .or. im /= 0.0) then
            t = h(j + 1, j) / slapy2(re, im)
            call cscal(kplusp - j + 1, conjg(t), h(j + 1, j), ncv)
            call cscal(min(j + 2, kplusp), t, h(1, j + 1), 1)
            call cscal(min(j + np + 1, kplusp), t, q(1, j + 1), 1)
            h(j + 1, j) = cmplx(real(h(j + 1, j)), 0.0)
        end if
    end do
    do i = 1, kev
        ! Final check for splitting and deflation. Use a standard test as in the QR algorithm
        ! REFERENCE: LAPACK subroutine zlahqr.
        ! Note: Since the subdiagonals of the compressed H are nonnegative real numbers, we take advantage of this.
        tst1 = cabs1(h(i, i)) + cabs1(h(i + 1, i + 1))
        if(tst1 == 0.0) tst1 = clanhs('1', kev, h, ncv, workl)
        if(real(h(i + 1, i)) <= max(ulp * tst1, smlnum)) h(i + 1, i) = 0.0
    end do
    ! Compute the (kev+1)-st column of (V*Q) and temporarily store the result in WORKD(N+1:2*N).
    ! This is needed in the residual update since we cannot GUARANTEE that the corresponding entry
    ! of H would be zero as in exact arithmetic.
    if (real(h(kev + 1, kev)) > 0.0) &
            call cgemv('N', n, kplusp, one, v, ldv, q(1, kev + 1), 1, zero, workd(n + 1), 1)
    ! Compute column 1 to kev of (V*Q) in backward order
    ! taking advantage of the upper Hessenberg structure of Q.
    do i = 1, kev
        call cgemv('N', n, kplusp - i + 1, one, v, ldv, q(1, kev - i + 1), 1, zero, workd, 1)
        call ccopy(n, workd, 1, v(1, kplusp - i + 1), 1)
    end do
    !  Move v(:,kplusp-kev+1:kplusp) into v(:,1:kev).
    call clacpy('A', n, kev, v(1, kplusp - kev + 1), ldv, v, ldv)
    ! Copy the (kev+1)-st column of (V*Q) in the appropriate place
    if (real(h(kev + 1, kev)) > 0.0) call ccopy(n, workd(n + 1), 1, v(1, kev + 1), 1)
    ! Update the residual vector:
    !    r <- sigmak*r + betak*v(:,kev+1)
    ! where
    !    sigmak = (e_{kev+p}'*Q)*e_{kev}
    !    betak = e_{kev+1}'*H*e_{kev}
    call cscal(n, q(kplusp, kev), resid, 1)
    if (real(h(kev + 1, kev)) > 0.0) call caxpy(n, h(kev + 1, kev), v(1, kev + 1), 1, resid, 1)
end
